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NUMERICAL INSTABILITY OF RESULTANT METHODS FOR 
MULTIDIMENSIONAL ROOTFINDING 

VANNI NOFERINI* AND ALEX TOWNSENDt 


Abstract. Hidden-variable resultant methods are a class of algorithms for solving multidimen¬ 
sional polynomial rootfinding problems. In two dimensions, when significant care is taken, they are 
competitive practical rootfinders. However, in higher dimensions they are known to miss zeros, cal¬ 
culate roots to low precision, and introduce spurious solutions. We show that the hidden variable 
resultant method based on the Cayley (Dixon or Bezout) matrix is inherently and spectacularly 
numerically unstable by a factor that grows exponentially with the dimension. We also show that 
the Sylvester matrix for solving bivariate polynomial systems can square the condition number of 
the problem. In other words, two popular hidden variable resultant methods are numerically unsta¬ 
ble, and this mathematically explains the difficulties that are frequently reported by practitioners. 
Regardless of how the constructed polynomial eigenvalue problem is solved, severe numerical dif¬ 
ficulties will be present. Along the way, we prove that the Cayley resultant is a generalization of 
Cramer’s rule for solving linear systems and generalize Clenshaw’s algorithm to an evaluation scheme 
for polynomials expressed in a degree-graded polynomial basis. 
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1. Introduction. Hidden variable resultant methods are a popular class of al¬ 
gorithms for global multidimensional rootfinding [niiziiiziiisiiMiiio]. They compute 
all the solutions to zero-dimensional polynomial systems of the form: 

. . .,Xd)\ 

: ^ ... ,Xd) G (1.1) 

\pdixi,...,Xd)) 

where d > 2 and pi,... ,pd are polynomials in xi,... ,Xd with complex coefficients. 
Mathematically, they are based on an elegant idea that converts the multidimensional 
rootfinding problem in (EH) into one or more eigenvalue problems |6]. At first these 
methods appear to be a practitioner’s dream as a difficult rootfinding problem is 
solved by the robust QR or QZ algorithm. Desirably, these methods have received 
considerable research attention from the scientific computing community [Tollllllaol 

H5] . 

Despite this significant interest, hidden variable resultant methods are notoriously 
difficult, if not impossible, to make numerically robust. Most naive implementations 
will introduce unwanted spurious solutions, compute roots inaccurately, and unpre- 
dictably miss zeros [5] . Spurious solutions can be removed by manually checking that 
all the solutions satisfy EH, inaccurate roots can usually be polished by Newton’s 
method, but entirely missing a zero is detrimental to a global rootfinding algorithm. 

The higher the polynomial degree n and the dimension d, the more pronounced 
the numerical difficulties become. Though our conditioning bounds do hold for small 
n and d, this paper deals with a worst-case analysis. Hence, our conclusions are 
not inconsistent with the observation that (at least when n and d are small) resultant 
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methods can work very well in practice for some problems. When d = 2 and real finite 
solutions are of interest, a careful combination of domain subdivision, regularization, 
and local refinement has been successfully used together with the Cayley resultant 
(also known as the Dixon or Bezout resultant) for large n [35]. This is the algorithm 
employed by Chebfun for bivariate global rootfinding [d^. Moreover, for d = 2, 
randomization techniques and the QZ algorithm have been combined fruitfully with 
the Macaulay resultant [27]. There are also many other ideas [11133]. However, these 
techniques seem to be less successful in higher dimensions. 

In this paper, we show that any plain vanilla hidden variable resultant method 
based on the Cayley or Sylvester matrix is a numerically unstable algorithm for solving 
a polynomial system. In particular, we show that the hidden variable resultant method 
based on the Cayley resultant matrix is numerically unstable for multidimensional 
rootfinding with a factor that grows exponentially with d. We show that for d = 2 
the Sylvester matrix leads to a hidden variable resultant method that can also square 
the conditioning of a root. 

We believe that this numerical instability has not been analyzed before because 
there are at least two other sources of numerical issues: (1) The hidden variable 
resultant method is usually employed with the monomial polynomial basis, which 
can be devastating in practice when n is large, and (2) Some rootfinding problems 
have inherently ill-conditioned zeros and hence, one does not always expect accurate 
solutions. Practitioners can sometimes overcome (1) W representing the polynomials 
Pi,... ,Pd in another degree-graded polynomial basiqj [H]- However, the numerically 
instability that we identify can be observed even when the roots are well-conditioned 
and for degree-graded polynomial basis (which includes the monomial, Chebyshev, 
and Legendre bases). 

We focus on the purely numerical, as opposed to symbolic, algorithm. We take 
the view that every arithmetic operation is performed in finite precision. There are 
many other rootfinders that either employ only symbolic manipulations [9] or some 
kind of symbolic-numerical hybrid m- Similar careful symbolic manipulations may 
be useful in overcoming the numerical instability that we identify. For example, it 
may be possible to somehow transform the polynomial system (ED into one that the 
resultant method treats in a numerical stable manner. 

This paper may be considered as a bearer of bad news. Yet, we take the opposite 
and more optimistic view. We are intrigued by the potential positive impact this 
paper could have on rootfinders based on resultants since once a numerical instability 
has been identified the community is much better placed to circumvent the issue. 

We use the following notation. The space of univariate polynomials with complex 
coefficients of degree at most n is denoted by C„ [a:], the space of d-variate polynomials 
of maximal degree n in the variables a:i,..., is denoted by Cn[xi, ..., x^], and if V 
is a vector space then the Cartesian product space V x • • • x V (d-times) is denoted 
by (V)'^. Finally, we use vec(V) to be the vectorization of the matrix or tensor V to 
a column vector (this is equivalent to V(:) in MATLAB). 

Our setup is as follows. First, we suppose that a degree-graded polynomial basis 
for C„[x], denoted by (j)o,... ,4>n, has been selected. All polynomials will be repre¬ 
sented using this basis. Second, a region of interest C is chosen such that 

where is the tensor-product domain D x • • • x D (d times), contains all the 
roots that would like to be computed accurately. The domain D C C can be a real 

polynomial basis {0o, ■ ■. for Cn[3;] is degree-graded if the degree of 0fe(a;) is exactly k 

for 0 < k < n. 
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interval or a bounded region in the complex plane. Throughout, we suppose that 
supjjgQ \4>k{x)\ = 1 for 0 < k < n, which is a very natural normalization. 

Our two main results are in Theorem EH and Theorem 14.61 Together they show 
that there exist pi,... ,pd in CH) such that 

Cond. no. of the eigenproblem Cond. no. of x* 

where R is either the Cayley (for any d > 2) or Sylvester (for d = 2) resultant matrix. 
Such a result shows that in the absolute sense the eigenvalue problem employed by 
these two resultant-based methods can be significantly more sensitive to perturbations 
than the corresponding root. Together with results about relative conditioning, we 
conclude that these rootfinders are numerically unstable (see Section Ell¬ 
in the next section we first introduce multidimensional resultants and describe 
hidden variable resultant methods for rootfinding. In Section El we show that the 
hidden variable resultant method based on the Cayley resultant suffers from numerical 
instability and in Section |4] we show that the Sylvester matrix has a similar instability 
for d = 2. In Section El we explain why our absolute conditioning analysis leads to 
an additional twist when considering relative conditioning. Finally, in Section El we 
present a brief outlook on future directions. 


2. Background material. This paper requires some knowledge of multidimen¬ 
sional rootfinding, hidden variable resultant methods, matrix polynomials, and con¬ 
ditioning analysis. In this section we briefly review this material. 


2.1. Global multidimensional rootfinding. Global rootfinding in high di¬ 
mensions can be a difficult and computationally expensive task. Here, we are con¬ 
cerned with the easiest situation where (EH has only simple finite roots. 

Definition 2.1 (Simple root). Let x* = (cc*, ...,G C‘^ be a solution to 
the zero-dimensional polynomial system EU. Then, we say that x* is a simple root 
of IHI]) if the Jacobian matrix J{x*) is invertible, where 


j(i>) = 


t(£") 




t-M"i 


dp±(^*\ 

dxi\^ 1. 


G C 


dxd 


( 2 . 1 ) 


If J{x*) is not invertible then the problem is ill-conditioned, and a numerically 
stable algorithm working in finite precision arithmetic may introduce a spurious solu¬ 
tion or may miss a non-simple root entirely. We will consider the roots of EB that 
are well-conditioned (see Proposition [2J]), finite, and simple. 

Our focus is on the accuracy of hidden variable resultant methods, not compu¬ 
tational speed. In general, one cannot expect to have a “fast” algorithm for global 
multidimensional rootfinding. This is because the zero-dimensional polynomial sys¬ 
tem in EB can potentially have a large number of solutions. To say exactly how 
many solutions there can be, we first must be more precise about what we mean by 
the degree of a polynomial in the multidimensional setting |38] . 

Definition 2.2 (Polynomial degree). A d-variate polynomial p{xi,... ,Xd) has 
total degree < n if 

d 

p{xi,...,Xd)= ^ 

iiH- hid<n k=l 
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for some tensor A. It is of total degree n if one of the terms ^ with ii + - ■ - + 1 ^ = 
n is nonzero. Moreover, p(xi,... ,Xd) has maximal degree < n if 

n d 

p{Xi,...,Xd)= ^ Aii,...,iciW<t>ik(.Xk) 

Zl,...,2d=0 k—1 

for some tensor A indexed by 0 < ii,... ,id < n. It is of maximal degree n if one of 
the terms with max(ii, ... ,id) = n is nonzero. 

Bezout’s Lemma says that if dm involves polynomials of total degree n, then 
there are at most solutions Chap. 3]. For polynomials of maximal degree we 
have the following analogous bound (see also O Thm. 5.1]). 

Lemma 2.3. The zero-dimensional polynomial system in dm, where pi,...,pd 
are of maximal degree n, can have at most dln‘^ solutions. 

Proof. This is the multihomogeneous Bezout bound, see [381 Thm. 8.5.2]. For 
polynomials of maximal degree n the bound is simply perm(n/d) = dln‘^, where Id is 
the d X d identity matrix and perm(A) is the permanent of A. □ 

We have selected maximal degree, rather than total degree, because maximal 
degree polynomials are more closely linked to tensor-product constructions and make 
later analysis in the multidimensional setting easier. We do not know how to repeat 
the same analysis when the polynomials are represented in a sparse basis set. 

Suppose that the polynomial system (ED contains polynomials of maximal degree 
n. Then, to verify that candidate points are solutions the polynomials pi,... ,pd 
must be evaluated, costing operations. Thus, the optimal worst-case complex¬ 
ity is For many applications global rootfinding is computationally unfeasible 

and instead local methods such as Newton’s method and homotopy continuation meth¬ 
ods |3] can be employed to compute a subset of the solutions. Despite the fact that 
global multidimensional rootfinding is a computationally intensive task, we still de¬ 
sire a numerically stable algorithm. A survey of numerical rootfinders is given in [441 
Chap. 5]. 

When d = 1, global numerical rootfinding can be done satisfactorily even with 
polynomial degrees in the thousands. Excellent numerical and stable rootfinders can 
be built using domain subdivision (3, eigenproblems with colleague or comrade ma¬ 
trices [23], and a careful treatment of dynamic range issues [7]. 

2.2. Hidden variable resultant methods. The first step of a hidden variable 
resultant method is to select a variable, say Xd, and regard the d-variate polynomials 
Pi,...,Pd in ED as polynomials in xi,..., Xd-i with complex coefficients that depend 
on Xd. That is, we “hide” Xd by rewriting Pk{xi ,..., Xd) for 1 < fc < d as 


n d—1 

Pk{xi,. . .,Xd-l,Xd) = Pk[Xd]ixi, . . .,Xd-l) = ^ 4>i.iXs), 

ii,...,id-i=0 s=l 

where {^o, • ■ •, </'«} is a degree-graded polynomial basis for C„[a:]. This new point of 
view rewrites ED as a system of d polynomials in d — 1 variables. We now seek all 
the S C such that Pi[x*d\,... ,Pd[x*d\ have a common root in Algebraically, 

this can be achieved by using a multidimensional resultant [20| Chap. 13]. 

Definition 2.4 (Multidimensional resultant). Let d > 2 and n > 0. A func¬ 
tional TZ : {Cn\xi,... ,Xd-i]Y C is a multidimensional resultant if, for any set 
of d polynomials qi,... ,qd € Cri[a;i,..., Xd-i], TZ-iqi, ■. ■, qu) is a polynomial in the 
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Fig. 1. Mathematically, the zeros of'R,{pi[X(i ],... ,Pd[^d]) component of the solutions 

to (II. ID . However, numerically the polynomial 'lZ{pi[xd], • • • ,Pdl^d]) numerically close to zero 

everywhere. Here, we depict the typical behavior of the polynomial 'lZ{pi[xd], ■ ■ • ,Pd[^d]) when d = 2, 
where the black dots are the exact zeros and the squares are the computed roots. In practice, it can 
be difficult to distinguish between spurious solutions and roots that are computed inaccurately. 

coejficients of qi,... ,qd and 'R-{qi ,..., qd) = 0 if and only if there exists an x* G 
such that qk{x*) = 0 for \ < k < d, where C denotes the extended complex 

Definition 12.41 defines 7^ up to a nonzero multiplicative constant [TTJ Thm. 1.6.1]. 
In the monomial basis it is standard to normalize TZ so that 7?.(x",... 1) = 

1 m Thm. 1.6.1(ii)]. For nonmonomial bases, we are not aware of any standard 
normalization. 

Assuming dni) only has finite solutions, if 7?. is a multidimensional resultant then 
for any xj £ C we have 

Ti{pi[xd],...,Pd[xd])=0 3(xJ,... ,x5_i) G s.t. pi(x*) = ---=pd(x*) = 0, 

where x* = (x*,... ,xj) G Thus, we can calculate the dth component of all the 
solutions of interest by computing the roots of TZ{pi[xd], ■ ■ ■ ,Pd[xd\) and discarding 
those outside of 17. In principle, since TZ{pi[xd], ■. ■ ,Pd[xd]) is a univariate polynomial 
in Xd it is an easy task. However, numerically, TZ is typically near-zero in large 
regions of C, and spurious solutions as well as missed zeros plague this approach in 
finite precision arithmetic (see Figure [1]). Thus, directly computing the roots of TZ 
is spectacularly numerically unstable for almost all n and d. This approach is rarely 
advocated in practice. 

Instead, one often considers an associated multidimensional resultant matrix 
whose determinant is equal to TZ. Working with matrices rather than determinants is 
beneficial for practical computations, especially when d = 2 [iTlESlEg]. Occasionally, 
this variation on hidden variable resultant methods is called numerically confirmed 
eliminants to highlight its improved numerical behavior [381 Sec. 6.2.2]. However, we 
will show that even after this significant improvement the hidden variable resultant 
methods based on the Cayley and Sylvester resultant matrices remain numerically 
unstable. 

Definition 2.5 (Multidimensional resultant matrix). Let d > 2, n > 0, TV > 1, 
and TZ be a multidimensional resultant (see Defintion\2.f^. A matrix-valued function 


^To make sense of solutions at infinity one can work with homogeneous polynomials m- 
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R : (C„[a;i,..., a^d-i])^ i^nxn ^ multidimensional resultant matrix associated 
with TZ if for any set of d polynomials qi, ... ,qd € C„[xi,..., x^-i] we have 

det {R{qi,...,qd)) = n{qi, . . . ,qd)- 


There are many types of resultant matrices including Cayley (see Section 131), 
Sylvester (see Section S]), Macaulay [17], and others [HI UHl 131] • In this paper we only 
consider two of the most popular choices: Cayley and Sylvester resultant matrices. 

Theoretically, we can calculate the dth component of the solutions by finding all 
the xj G C such that det(i?(pi[xj],... ,p£;[x5])) = 0. In practice, our analysis will 
show that this dth component cannot always be accurately computed. 

Each entry of the matrix R(pi[xd\, ■ ■ ■ ,Pd[xd\) is a polynomial in Xd of finite 
degree. In linear algebra such objects are called matrix polynomials (or polynomial 
matrices) and finding the solutions of det(i?(pi[x£;],... ,Pd[xd\)) = 0 is related to a 
polynomial eigenproblem [31311113]. 

2.3. Matrix polynomials. Since multidimensional resultant matrices are ma¬ 
trices with univariate polynomial entries, matrix polynomials play an important role 
in the hidden variable resultant method. A classical reference on matrix polynomials 
is the book by Gohberg, Lancaster, and Rodman [22] . 

Definition 2.6 (Matrix polynomial). Let N >1 and K >0. We say that P{X) 
is a (square) matrix polynomial of size N and degree K if P{X) is an N x N matrix 
whose entries are univariate polynomials in X of degree < K, where at least one entry 
is of degree exactly K. 

In fact, since is a zero-dimensional polynomial system it can only have a 

finite number of isolated solutions and hence, the matrix polynomials we consider are 
regular [22] . 

Definition 2.7 (Regular matrix polynomial). We say that a square matrix 
polynomial P{X) is regular if det{P{X)) ^ 0 for some A G C. 

A matrix polynomial P{X) of size N and degree K can be expressed in a degree- 
graded polynomial basis as 


K 

P{X)=J2A^M^), A, gC^^^. (2.2) 

i=0 

When the leading coefficient matrix Ak in (12.2|) is invertible the eigenvalues of P{X) 
are all finite, and they satisfy det(P(A)) = 0. 

Definition 2.8 (Eigenvector of a regular matrix polynomial). Let P{X) be a 
regular matrix polynomial of size N and degree K. If X G C is finite and there exists 
a non-zero vector v G such that P{X)v = 0 (resp. u^P(A) = 0), then we say 

that V is a right (resp. left) eigenvector of P{X) corresponding to the eigenvalue A. 

For a regular matrix polynomial P{X) we have the following relationship between 
its eigenvectors and determinant [22]: For any finite A G C, 

det(P(A))=0 ^ 3 i;G C^^^\{0}, P{X)v = 0. 

In multidimensional rootfinding, one sets P{X) = i?(pi[A],... ,pd[A]) and solves 
det(P(A)) = 0 via the polynomial eigenvalue problem P(X)v = 0. There are various 
algorithms for solving P{X)v = 0 including linearization [H] [311 133], the Ehrlich- 
Aberth method [3 mi mi, and contour integration [3]. However, regardless of how 
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the polynomial eigenvalue problem is solved in finite precision, the hidden variable 
resultant method based on the Cayley or the Sylvester matrix is numerically unstable. 

For the popular resultant matrices, such as Cayley and Sylvester, the first d — \ 
components of the solutions can be determined from the left or right eigenvectors of 
R{pi[x*^,... ,pd[x’^]). For instance, if linearization is employed, the multidimensional 
rootfinding problem is converted into one (typically very large) eigenproblem, which 
can be solved by the QR or QZ algorithm. Practitioners often find that the computed 
eigenvectors are not accurate enough to adequately determine the d — 1 components. 
However, the blame for the observed numerical instability is not only on the eigen¬ 
vectors, but also the eigenvalues. Our analysis will show that the dth component may 
not be computed accurately either. 

2.4. Conditioning analysis. Not even a numerically stable algorithm can be 
expected to accurately compute a simple root of (HID if that root is itself sensitive 
to small perturbations. Finite precision arithmetic almost always introduces roundoff 
errors and if these can cause large perturbations in a root then that solution is ill- 
conditioned. 

The absolute condition number of a simple root measures how sensitive the loca¬ 
tion of the root is to small perturbations in pi,... ,pd. 

Proposition 2.9 (The absolute condition number of a simple root). Let x* = 
{x\,... ,x*^) G be a simple root of (HID. The absolute condition number of x* 
assoeiated with rootfinding is ||T(a:*)“^|| 2 , i-e., the matrix 2-norm of the inverse of 
the Jacobian. 

Proof. See [35]. □ 

As a rule of thumb, a numerically stable rootfinder should be able to compute a 
simple root x* G to an accuracy of 0(max(|| J(x*)“^|| 2 , l)u), where u is the unit 
machine roundoff. In contrast, regardless of the condition number of s*, a numerically 
unstable rootfinder may not compute it accurately. Worse still, it may miss solutions 
with detrimental consequences. 

A hidden variable resultant method computes the dth component of the solutions 
by solving the polynomial eigenvalue problem R{pi[xd], ■. ■ ,Pd[xd])v = 0. The follow¬ 
ing condition number tells us how sensitive an eigenvalue is to small perturbations in 
R [3S1 (12)] (also see [33]1: 

Definition 2.10 (The absolute condition number of an eigenvalue of a regular 
matrix polynomial). Letx*^ G C be a finite eigenvalue of R{xd) = R(jii[xd], ■. ■ ,pd[xd\)- 
The condition number of x’^ associated with the eigenvalue problem R{xd)v = 0 is 

K{Xd, R) = lim sup I - min \xd — x'^\ : det(R(xd)) = 01 , (2.3) 

£->0+ (e J 

where the supremum is taken over the set of matrix polynomials R(xd) such that 
max,rdGn \\Rixd) - R{xd)\\2 < £• 

A numerical polynomial eigensolver can only be expected to compute the eigen¬ 
value ccj satisfying R(x*^)v = 0 to an accuracy of 0(max(K(a;J, i?), l)u), where u is 
unit machine roundoff. We will be interested in how k{x’^,R) relates to the condition 
number, || J(a:*)“^ || 2 , of the corresponding root. 

It can be quite difficult to calculate k{x*^,R) directly from p.3ll . and is usually 
more convenient to use the formula below. (Related formulas can be found in [351 
Thm. 1] for symmetric matrix polynomials and in [421 Thm. 5] for general matrix 
polynomials.) 


Lemma 2.11. Let R{xd) be a regular matrix polynomial with finite simple eigen¬ 
values. Let € C he an eigenvalue of R{xd) with corresponding right and left eigen¬ 
vectors v,w G Then, we have 


K{x*d,R) 


Ikihllwlh 

\w'^R'{xd)v \’ 


where R'(xd) denotes the derivative of R with respect to Xd- 

Proof. The first part of the proof follows the analysis in [42]. Let R{xd) be a 
regular matrix polynomial with a simple eigenvalue ccj G C and corresponding right 
and left eigenvectors v,w G A perturbed matrix polynomial R{x) = R{x) + 

AR{x) will have a perturbed eigenvalue Xd and a perturbed eigenvector v = v -\- 6v 
such that R{xd)v + AR{xd)v = 0, where || Ai?(a ;)||2 < £■ 

Expanding, keeping only the first order terms, and using R(x*^v = 0 we obtain 

{xd - x*d)R'{x*d)v + R{x*d)5v + AR{x*d)v = 0{e^). 


Multiplying by on the left, rearranging, and keeping the first order terms, we 
obtain 


, w'^AR{x*)v 
= w^R'ix^v ’ 

where the derivative in R' {x*^) is taken with respect to Xd. Thus, from (12.3p we see 
that 


K{x*d,R) < 


Ikll2||w||2 

\w'^R'{x*^)v\' 


(2.4) 


We now show that the upper bound in (12.41) can be attained. Take AR{xd) = 
ewu'^/(||u|| 2 ||w|| 2 )- Then, max,,^ga || Ai?(a:d )||2 = e and 


w'^AR{x*d)v ^ lkl|2||w||2 

w'^R'{x*j)v w'^ R'{x*jfv' 


The result follows by Definition 12.101 □ 

For the Cayley resultant matrix (see Section |3|), we will show that K 2 {x*j^,R) can 
be as large as ||J (^*)“^||2 (see Theorem 13.71) . Thus, there can be an exponential 
increase in the conditioning that seems inherent to the methodology of the hidden 
variable resultant method based on the Cayley resultant matrix. In particular, once 
the polynomial eigenvalue problem has been constructed, a backward stable numerical 
eigensolver may not compute accurate solutions to (HD. 

We now must tackle the significant challenge of showing that the Cayley and 
Sylvester resultant matrices do lead to numerical unstable hidden variable resultant 
methods, i.e., for certain solutions x* the quantity K 2 (x^, R) can be much larger than 

3. The Cayley resultant is numerically unstable for multidimensional 
rootfinding. The hidden variable resultant method when based on the Cayley re¬ 
sultant [12] finds the solutions to (HD by solving the polynomial eigenvalue problem 
given by Rcayiey{xd)v = 0, where Rcayiey{xd) is a certain matrix polynomial. To 
define it we follow the exposition in m and first introduce a related Cayley function 

f Cayley ■ 
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Definition 3.1 (Cayley function). The Cayley function associated with the poly¬ 
nomials qi, ■ ■ ■ ,qd ^ ..., Xd-i] is a multivariate polynomial in 2d —2 variables, 

denoted by fcayiey = fCayley(qi,- ■ ■, qd), and is given by 


fCayley — d^t 


fqiisi, S2, ■ ■ ■ , Sd-l) ... qdisi, S2, ■ ■ ■ , Sd-l)\ 
qiitl, S2, ■ ■ ■ , Sd-l) ... qd(ti, S2, ■ ■ ■ , Sd-l) 

^qi{ti,t2, ■ ■ ■ ,td-l) ... qd{ti,t2, ■ . ■ ,td-l) j 


(3-1) 


In two dimensions the Cayley function (also known as the Bezoutian function 
takes the more familiar form of 


fCayley — 


Si - tl 


■ det 


gi(si) g 2 (si)' 

,qi{ti) q2{ti), 


qi{si)q 2 {ti) - q 2 {si)qi{ti) 
Si - tl 


which is of degree at most n — 1 in si and ti. By carefully applying Laplace’s formula 
for the matrix determinant in (13.11) , one can see that fcayiey is a polynomial of degree 
Tfc < fcn — 1 in Sk and td-k for 1 < fc < d — 1. 

Note that fcayiey is not the multidimensional resultant (except when = 0 for 
all k). Instead, fcayiey is a function that is a convenient way to define the Cayley 
resultant matrix. 

Let {(j)o,4>i, ■ ■ ■,} be the selected degree-graded polynomial basis. The Cayley 
resultant matrix depends on the polynomial basis and is related to the expansion 
coefficients of fcayiey in a tensor-product basis of {(fo, (pi ,...That is, let 


Tl Td-1 Td-1 n d—1 d—1 

fCay ley = ^ ^ ^ ^ ^ ^ ^*1 . .,id-i ,ti,.. .,id-l ('Sfc ) ) (3-2) 

n=0 id_i=0ji=0 jd-l=0 k=l k=l 

be the tensor-product expansion of the polynomial fcayiey, where A is a tensor of 
expansion coefficients of size (ri -|- 1) x • • • x (r^-i -(- 1) x (r^-i -I- 1) x • • • x (ri -|- 1). 
The Cayley resultant matrix is the following unfolding (or matricization) of A [361 
Sec. 2.3]: 

Definition 3.2 (Cayley resultant matrix). The Cayley resultant matrix asso¬ 
ciated with qi,...,qd € Cn[xi, ... ,Xd-i] with respect to the basis {(po, pi,...,} is 
denoted by Rcayiey and is the (nfe=i('^fc + 1 )) ^ (nfc=i(''‘fc + 1 )) matrix formed by 
the unfolding of the tensor A in (13.21) . This unfolding is often denoted by Arxr, where 
r ={!,..., d-1} and c = {d,...,2d-2} [31 Sec. 2.3]. 

For example, when Tk = kn—1 for l<A:<d—Iwe have for 0 < ik,jd-k < kn—1 

C d-l d-l 

'y i)hfcn 1 y ^ jd—k 

k—l k—1 

This is equivalent to N = factorial (d-l) (d-l) ; R = reshape(A, N, N) ; in 

MATLAB, except here the indexing of the matrix Rcayiey starts at 0. 

For rootfinding, we set qi = pi[xd], ■ ■ ■ ,qd = Pd[xd\ (thinking of Xd as the “hid¬ 
den” variable). Then, Rcayiey = Rcayiey(xd) is a square matrix polynomial (see 
Section HH). If all the polynomials are of maximal degree n, then Rcayiey is of size 


id^k-, 

{d-k)\ 


= A. 
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{d — l)\n‘^ ^ and of degree at most dn. The fact that (d — l)!n'^ ^ x dn = d\nf' is 
the maximum number of possible solutions that (HH) can possess (see Lemma lOl) is 
a consequence of Rcayiey being a resultant matrix. In particular, the eigenvalues of 
Rcayieyixd) are the dth components of the solutions to (11.11) and the remaining d—\ 
components of the solutions can in principle be obtained from the eigenvectors. 

It turns out that evaluating fcayiey at is equivalent to a matrix-vector 

product with Rcayiey ■ This relationship between Rcayiey and fcayiey will be essential 
in Section for understanding the eigenvectors of Rcayiey 

Lemma 3.3. Let d > 2, t* G and fcayiey o,nd Rcayiey be the Cayley 

function and matrix associated with qi,... ,qd G ..., Xd-i], respectively. If V 

is the tensor satisfying 'Pjki^k) ^ — jd-k < Tfe, then we have 

Rcayleyyec{v) = Vec(y), 

where Y is the tensor that satisfies 


Ti Td-l d—1 

fCayleyjsi, . . . , Sd-1, t^, . . . , t^_i) = ,■ ■ ■ ,i<i-1 ) • 

ii—0 id-1—0 k—1 


Proof. The matrix-vector product Rcayiey'veclV) = vec{Y) is equivalent to the 
following sums: 


Ti d—1 

jd-l=0 k=l 

for some tensor Y. The result follows from (13.21) . □ 

3.1. The Cayley resultant as a generalization of Cramer’s rule. In this 
section we show that for systems of linear polynomials, i.e., of total degree 1, the 
Cayley resultant is precisely Cramer’s rule. We believe this connection is folklore, but 
we have been unable to find an existing reference that provides a rigorous justification. 
It gives a first hint that the hidden variable resultant method in full generality may 
be numerically unstable. 

Theorem 3.4. Let A be a matrix of size d x d, x = (a;i,..., Xd)’^, and b a vector 
of size d X 1. Then, solving the linear polynomial system Ax + b = 0 by the hidden 
variable resultant method based on the Cayley resultant is equivalent to Cramer’s rule 
for calculating Xd. 

Proof. Let Ad be the last column of A and B = A—AdC^ Pbe^, where Cd is the dth 
canonical vector. Recall that Cramer’s rule computes the entry Xd in Ax = —b via the 
formula Xd = — det(i?)/det(A). We will show that for the linear polynomial system 
Ax -I- 6 = 0 we have fcayiey = det{B) + Xd det(A). Observe that this, in particular, 
implies that (since fcayley has degree 0 in S^, U for all i) fcayley = RCayley = R-Cayley 
Hence, the equivalence between Cramer’s rule and rootfinding based on the Cayley 
resultant. 

First, using (ED, we write fcayiey = det(M)/det(C) where the matrices M and 


E- 

ii=o 


11 


V are 


51 ti ti ... ti 

52 32 t2 ... t2 


v = 


M = BV + XdAdC^, 


Sd-l 

1 


Sd-l 

1 


Sd-l 

1 


id-l 

1 


where e is the d x 1 vector of all ones. (It can be shown by induction on d that 
det(y)=nti(s-i*). as required.) Using the matrix determinant lemma, we have 

det(M) = det(i?) det(U) + Xde^ &d]{BV)Ad, 


where adj(i3U) is the algebraic adjugate matrix of BV. Now, recall that adj(i3U) = 
adj(U) adj(i3) and observe that adj(U) = det(U)eJ. Hence, we obtain 


det(M) 

det(U) 


det(H) + Xd{e^ adj(H)Hd). 


Using ej adj(i3)6 = det(i?) and the matrix determinant lemma one more time, we 
conclude that 


det(H) = det(H) + ej a.d]{B)Ad - ej adj(H)6 = ej &d]{B)Ad. 

Thus, fCayley = det(i?) + Xddet{A) and the resultant method calculates Xd via 
Cramer’s formula. □ 

It is well-known in the literature that Cramer’s rule is a numerically unstable 
algorithm for solving Ax = b [531 Sec. 1.10.1]. Thus, Theorem 13.41 casts significant 
suspicion on the numerical properties of the hidden variable resultant method based 
on the Cayley resultant. 

3.2. The eigenvector structure of the Cayley resultant matrix. Ulti¬ 
mately, we wish to use Lemma 12.111 to estimate the condition number of the eigen¬ 
values of the Cayley resultant matrix. To do this we need to know the left and right 
eigenvectors of Rcayiey The following lemma shows that the eigenvectors of Rcayiey 
are in Vandermonde fornj^. To show this we exploit the convenient relationship be¬ 
tween evaluation of fcayiey and matrix-vector products with Rcayiey 

Lemma 3.5. Suppose that x* = (a:^,... ,x^) G is a simple root of (11.111 . Let 
V and W be tensors of size (jd-i -|- 1) x • • • x (ti -I- 1) and (ti -|- 1) x • • • x {xd-i + 1), 
respectively, defined by 

d-l 

^1, -Jd-l = 0 < Jfc < Td-k 

fe=i 

and 

d-l 

0 < ifc < Tk- 

k=l 


®In one dimension we say that an iV x 1 vector v is in Vandermonde form if there is an x € C 
such that Vi = for 0 < i < Ai — 1. In higher dimensions, the vector vec(yl) is in Vandermonde 

form if = nf=i (^fc) for some xi,... ,Xd G C. 
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Then, the vectors vec(l^) and vec{W) are the right and left eigenvectors of the matrix 
Rc ayiey(pi[x*,j\, ■ ■ ■ ^Pd{x*^) that Correspond to the eigenvalue xj. 

Proof Let fcayiey = fCayleyi.Piix'fl, ■. ■,Pd[x^]) be the Cayley function associated 
withpi[a;5],... ,Pd[x*j]. From dMl we find that fcayiey{si, ..., Sd-i,xl ,... ,a:2_ J = 0 
because the determinant of a matrix with a vanishing last row is zero. Moreover, by 
Lemma 13.31 we have 

Ti Td-1 d—1 

0 = fcayleyi^l^ • ■ - ; . . . , ^ ^ ^ ^ 4^ik i^k)- 

ii—0 k—1 


Since {(fo, 4>i, ■ ■ ■,} is a polynomial basis we must conclude that F = 0, and hence, 
Rcayiey{x*^)v = 0 with V = vec(F). In other words, n is a right eigenvector of Rcayiey 
corresponding to the eigenvalue (see Definition 12.8p . 

An analogous derivation shows that vec(IF) is a left eigenvector of Rcayiey □ 

3.3. On the generalized Rayleigh quotient of the Cayley resultant ma¬ 
trix. To bound k{x’^, Rcayiey) we need to bound the absolute value of the generalized 
Rayleigh quotient of R'cayieyixd) (see Lemma [2. Ill) , whenever x* G is such that x’^ 
is a simple eigenvalue of Rcayiey{xd), i-e., there are no other solutions to (11.11) with 
the same dth component. In a similar style to the proof of Lemma 13.51 we show this by 
exploiting the relation between evaluating the derivative of fcayiey and matrix-vector 
products with R'cayleyi^d)- 

Theorem 3.6. Let pi,... ,pd be the polynomials in (HD, X* G a solution 
of dHD, and fcayiey{xd) the Cayley function associated with qi = pi[xd], ■ ■ ■ ,qd = 
Pd[xd\- We have 


f'cavlev{x*d) 


l<k<d-l 


det(J(4)), 


where J{x*) is the Jacobian matrix in (12.111 . That is, f'cayieyixd) evaluated at Sk = 
tk = for 1 < k < d — 1 is equal to the determinant of the Jacobian. 

Proof. Recall from that fcayiey{xd) is a polynomial in si,...,Sd_i and 

ti, ..., td-i written in terms of a matrix determinant, and set qi = pi[xd\,.. ■ ,qd = 
Pd[xd]- The determinant in (13.1|) for fcayiey{xd) can be expanded to obtain 


fCayley{,Xd) — ^ ^ ( 1 ) Pai\Xd\{f'l’: ■ • ■ • ■ • , 

lli=l (TGSd i=l 

where Sd is the symmetric group of {l,...,(i} and (—1)“^ is the signature of the 
permutation a. When we evaluate fcayiey{xd) aJ Sk = tk = x"^ ior 1 < k < d — \ 
the denominator vanishes, and hence, so does the numerator because fcayiey{xd) is 
a polynomial. Thus, by L’Hospital’s rule, fcayieyixd) evaluated Sk = tk = x^. for 
l<A:<d—lis equal to 


dsi 


dd 

■■dsd-idxd 


d 

'y ^ ( 1) Pat [2^(i](^l) ■ • ■ J ti—i, Si, ... , Sd—l) 

aCiSd i=l 


(3.3) 


evaluated at Sk = x\, tk = x^., and Xd = x*y In principle, one could now apply the 
product rule and evaluate the combinatorially many terms in (13.31) . Instead, we note 
that after applying the product rule a term is zero if it contains Pa^ (x*) for any a G Sd 
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and 1 <i < d (since a:* is a solution to (ED)- There are precisely d partial derivatives 
and d terms in each product so that any nonzero term when expanding 13.31 has each 
Pk differentiated precisely once. Finally, note that for each 1 < fc < d — 1 only the 
\ <i <k terms in the product depend on Sk- Hence, from (13.31) we obtain 


fcayleyi^d) 


cGSd i=l 


l<fe<d-l 


dpa 

dx, 


-fe*)- 


The result follows because the last expression is the determinant of the Jacobian 
matrix evaluated at x*. □ 

As a consequence of Theorem 13.61 we have the following unavoidable conclusion 
that mathematically explains the numerical difficulties that practitioners have been 
experiencing with hidden variable resultant methods based on the Cayley resultant. 

Theorem 3.7. Let d> 2. Then, there exist pi, ■. ■ ,pd in EU with a simple root 
X* S such that 


k{Xj, RCayley) > ) II2 


and ||J(ai*) ^||2 > 1- Thus, an eigenvalue of Rcayiey{xd) can be more sensitive to 
perturbations than the corresponding root by a factor that grows exponentially with d. 


Proof Using Lemma 13.31 Theorem 13.61 has the following equivalent matrix form: 

R'cayley{x*d)v = det(J(x*)), 

where v = vec(U), w = vec(lU), and V and W are given in Lemma 15751 Since (fo = 1, 
we know that ||'(;||2 > 1 and ||rc ||2 > 1- Hence, by Lemma [2.Ill 

nix^, Rcayley) > | det( J(x*)) 

Denoting the singular values [IHl Sec. 7.3] of the matrix J{x*) by ai , select pi,...,pd 
and X* G such that |det(J(x*))| = 0^=1= ^d- Such polynomial systems do 
exist, for example, linear polynomial systems where Mx — Mx* = 0 and M is a matrix 
with singular values (Ti = <72 = ■•■ = crd. To ensure that ||J (®*)“^||2 > 1 we also 
require < 1. Then, we have 

d 

K{Xd, Rcayley)~^ < |det(J(x*))| 

i=l 


The result follows. □ 

Example 3.8. Let Q be a dx d orthogonal matrix, QQ^ = Id, having elements 
Qij for i, j = 1,..., d, and let a < 1. Consider the system of polynomial equations 

d 

Pi = xf + cr'^qijXj =0, i = l,...,d. 
i=i 

The origin, x* = 0 £ C'^, is a simple root of this system of equations. The Jacobian 
of the system at 0 is J = aQ, and hence, the absolute conditioning of the problem is 
||J“^|| = . Constructing the Cayley resultant matrix polynomial in the monomial 

basis, one readily sees that for this example the right and left eigenvectors for the 
eigenvalue x’^ = 0 satisfy ||u|| = ||■u;|| = 1. As a consequence, k{x’^, Rcayiey) = cr~‘^. 
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We emphasize that this numerical instability is truly spectacular, affects the ac¬ 
curacy of a;J, and can grow exponentially with the dimension d. 

Moreover, Theorem 1X71 holds for any degree-graded polynomial basis selected 
to represent pi,...,pd as long as (/jq = 1. In particular, the associated numerical 
instability cannot be resolved in general by a special choice of polynomial basis. 

Theorem 1X71 is pessimistic and importantly does not imply that the resultant 
method always loses accuracy, just that it might. In general, one must know the 
solutions to CH) and the singular values of the Jacobian matrix to be able to predict 
if and when the resultant method will be accurate. 

One should note that Theorem 13.71 concerns absolute conditioning and one may 
may wonder if a similar phenomenon also occurs in the relative sense. In Section [5] 
we show that the relative conditioning can also be increased by an exponential factor 
with d. 

4. The Sylvester matrix is numerically unstable for bivariate rootfind¬ 
ing. A popular alternative in two dimensions to the Cayley resultant matrix is the 
Sylvester matrix m Chap. 3], denoted here by Rsyiv We now set out to show that 
the hidden variable resultant based on Rsyiv is also numerically unstable. However, 
since d = 2 the instability has only a moderate impact in practice as the conditioning 
can only be at most squared. With care, practical bivariate rootfinders can be based 
on the Sylvester resultant [31] though there is the possibility that a handful digits are 
lost. 

A neat way to define the Sylvester matrix that accommodates nonmonomial poly¬ 
nomial bases is to define the matrix one row at a time. 

Definition 4.1 (Sylvester matrix). Let qi and q 2 be two univariate polynomials 
in^2n[xi] of degree exactly ti andT 2 , respectively. Then, the Sylvester matrix Rsyiv G 
([;;(Ti+T 2 )x(ri+r 2 ) associated with qi and q 2 is defined row-by-row as 

Rsyiv (*,:)= 0 < i < r2 - 1, 

where is the row vector of coefficients such that qi{x)(j)i{x) = ^k"^4‘k{x) 

and 


Rsyiv{'i + r 2 , ■■) 0 < Z < Ti - I, 


where is the row vector of coefficients such that q 2 {x)(pi{x) = ^ '^k‘^4‘k{x). 


In the monomial basis, i.e., 4>k{x) = x^, Definition 14.11 gives the Sylvestei0 matrix 
of size (ti -I- r 2 ) x (n -|- T 2 ) as [T5| Chap. 3]: 


Rsylv 


/ao 

Oi 


^Tl 


\ 



ao 

ai . 

Qjti 

\ 

bo 

bi 


^T2 


\ 

V 


bo 

h - 

br2 j 



T 2 rows 


Ti rows 


(4.1) 


“^Variants of (14.11 1 include its transpose or a permutation of its rows and/or columns. Our 
analysis still applies after these aesthetic modifications with an appropriate change of indices. We 
have selected this variant for the convenience of indexing notation. 
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where qi{x) = Y2 =q o-kx'" and q 2 {x) = Y.'k^a ^kx’"- 

4.1. A generalization of Clenshaw’s algorithm for degree-graded poly¬ 
nomial bases. Our goal is to use Lemma [2 .111 to bound the condition number of the 
eigenvalues of the Sylvester matrix. It turns out the right eigenvectors of Rsyiv are 
in Vandermonde form. However, the left eigenvectors have a more peculiar structure 
and are related to the byproducts of a generalized Clenshaw’s algorithm for degree- 
graded polynomial bases (see Lemma [4.411 . We develop a Clenshaw’s algorithm for 
degree-graded bases in this section with derivations of its properties in Appendix 

The selected polynomial basis (j)o, (j)i,..., is degree-graded and hence, satisfies a 
recurrence relation of the form 

k 

4)k+i{x) = (akX + Pk)(l>k{x)+'^'yk,j(l>j-i{x), k>l, (4.2) 

i=i 

where 4>i{x) = (agx + /3o)(/>o(a^) and (fio(x) = 1. If t/fo, , is an orthogonal poly¬ 

nomial basis, then (14.21) is a three-term recurrence and it is standard to employ Clen¬ 
shaw’s algorithm m to evaluate polynomials expressed as p{x) = Y^2=o^k'Pk{,x). 
This procedure can be extended to any degree-graded polynomial basis. 

Let p{x) be expressed as p{x) = '^k4>kix), where (j)o, ■ ■ ■,'Pn is a degree- 

graded polynomial basis. One can evaluate p{x) via the following procedure: Let 
bn+i[p]{x) = 0, and calculate bn[p]{x ),..., bi[p]{x) from the following recurrence rela¬ 
tion: 

n—1 

bk[p]ix) = Qk + {akX + Pk)bk+i[p]{x) + ^ Jj^k+ibj+i[p]{x), 1 < k < n. (4.3) 

j=k+l 

We refer to the quantities bi[p\(x ),..., bn+i[p]{x) as Clenshaw shifts (in the monomial 
case they are called Horner shifts [16]). The value p{x) can be written in terms of the 
Clenshaw shiftfl 

Lemma 4.2. Let n be a positive integer, x G C, (po,... a degree-graded 
basis satisfying (14.2p . p{x) = J2k=o ^kpkix), and bn+i[p]{x ),..., bi[p\{x) the Clenshaw 
shifts satisfying (14.31) . Then, 

n—1 

p{x) = aocpoix) -G (l)iix)bi[p]{x) + '^ji^ibi+i[p]{x). (4.4) 

2 = 1 


Proof See Appendix □ 

Clenshaw’s algorithm for degree-graded polynomial bases is summarized in Fig¬ 
ure [2| We note that because of the full recurrence in (14.31) the algorithm requires 
0{n^) operations to evaluate p{x). Though this algorithm may not be of significant 
practical importance, it is of theoretical interest for the conditioning analysis of some 
linearizations from the so-called Li- or L 2 -spaces [3l| when degree-graded bases are 
employed |5i] . 

There is a remarkable and interesting connection between Clenshaw shifts and the 
quotient {p(x) — p{y))/{x — y), which will be useful when deriving the left eigenvectors 
of Rsylv ■ 

®Note that, although Lemma 14.21 is stated in a general form and holds for any degree-graded 
basis, in this paper we fix the normalization maxj;gn \(pj{x)\ = 1, that implies in particular 0o = 1 
simplifying 114. 211 . 
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Clenshaw’s algorithm for degree-graded polynomial bases 

Let 00 , 01 ,, satisfy and p{x) = J2k=o o-k4>k{x). 

Set hn+i[p]{x) = 0. 

for A: = n,n — 1,...,1 do 

hk[p]{x) = ak + {akX +l3k)hk+i[p\{.x) + YTjZk+ilj,k+ibj+i[p\{x) 

end 

p{x) = 0000 (a;) + (l>i{x)bi[p]{x) + YJjZl Ijpbj+iivM^)- 

Fig. 2. Clenshaw’s algorithm for evaluating polynomials expressed in a degree-graded basis. 


Theorem 4.3. 


With the same set up as Lemma \4^ we have 

^ a^b^+l[p]{y)(|)^{x), x^y 

nr. — 7/ 


and 


n— 1 

= ^azb,+i[p]{x)(l)i(x). 
2=0 


(4.5) 


(4.6) 


Proof. See Appendix □ 

The relation between the derivative and Clenshaw shifts in (14.61) has been noted 
by Skrzipek for orthogonal polynomial bases in m, where it was used to construct a 
so-called extended Clenshaw’s algorithm for evaluating polynomial derivatives. Using 
Theorem 14.31 and [S^ an extended Clenshaw’s algorithm for polynomials expressed in 
a degree-graded basis is immediate. 

4.2. The eigenvector structure of the Sylvester matrix. We now set qi = 

Pi[x 2 ] and q 2 = _P 2 [a: 2 ] (considering X 2 as the hidden variable), and we are interested 
in the eigenvectors of the matrix polynomial Rsyivix^), when (xi,X 2 ) is a solution 
to (|l.ll) when d = 2. It turns out that the right eigenvectors of Rsyivix^) are in 
Vandermonde form, while the left eigenvectors are related to the Clenshaw shifts (see 
Section HD) . 

Lemma 4.4. Suppose that x* = (ccj,^^) is a simple root of dm and that 
Pi[x 2 ] and P 2 [x 2 ] are of degree ti and T 2 , respectively, in xi. The right eigenvec¬ 
tor of Rsyiv{x 2 ) corresponding to the eigenvalue is 

Vk = (l>kix*i), 0 < fc < Ti -I- T2 - 1, 


and the left eigenvector is defined as 


Wi = 


-aibi+i[q2]{x\), 

at-r^b,_r^+i[qi]{xl), 


0 < i < r2 — 1, 

T2 < i < Ti T2 — I, 


where qj = Pj[x 2 ] and bk[qj]{x’l) are the Clenshaw shifts with respect to {0o, 0i,...,}, 
while the coefficients Ui are defined as in (I4.2|l . 

Proof. By construction we have, for 0 < i < r 2 — 1, 

T1+T2 — 1 

RSylv{i,-)v= X! (^ 2 )(l^k{xl) = qi{xl)(j)i{xl) = 0 
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and, for 0 < i < Ti — 1, 


T1+T2 — 1 

Rsylv {i + T 2 ,:)v= ^ Yl’‘^{x*2)4>k{x\) = q2{xl)4>i{xl) = 0. 

k=0 

Thus, u is a right eigenvector of Rgyivix^) corresponding to the eigenvalue x\. 

For the left eigenvector, first note that for any vector $ of the form = 4)k{x) 
for 0 < fc < Ti + r 2 — 1 we have by Theorem 14.31 


T2 — 1 Ti — 1 

■uFRsyivixl)^ = - ^ aihi+i[q2]{x\)(j)i{x)qi{x) + ^ aibi+i[qi]{xl) 4 >i{x)q 2 {x) 
2—0 2—0 

q2{x)-q2{xl) , , , qi{x)-qi{xl) 

- -qi{x)-{ -;- q2{x) 


X — X 
<12{X) 


qi{x) + ^^^^\ q 2 {x) = 0, 


where the second from last equality follows because qi{x\) = q 2 {x\) = 0. Since (j4.3p 
holds for any x and {i^O; ; </^’ti+t 2 -i} is a basis of CT-i+r 2 -iN: we deduce that 

iFRsyivixF) = 0, and hence, w is a left eigenvector of Rsyiv corresponding to the 
eigenvalue □ 

4.3. On the generalized Rayleigh quotient of the Sylvester matrix. To 

bound K{Rsyiv,x’^) we look at the absolute value of the generalized Rayleigh quo¬ 
tient of R'syivix 2 ), whenever x* is such that xj is a simple eigenvalue of Rsyiv{x 2 )- 
Lemma [4.4l allows us to show how the generalized Rayleigh quotient of R'syivFV) 
relates to the determinant of the Jacobian. 

Lemma 4.5. With the same assumptions as in Lemma \4-4\ we have 

\w^R'syivF^Fl ^ |det(J(x*))| 

lkl |2||«^||2 - Ik ||2 ’ 


where w and v are the left and right eigenvectors of Rsyiv, respectively, and J(x*) 
the Jacobian matrix in m- 

Proof. By Lemma |4.4l we know the structure of v and w. Hence, we have 


R'syivFVtv = -^a,h+i[q2\{xl)(l)i{x\)^^{x\) p'^a,bi+i[qi]{x\)(t),{x\)^{xl) 


2=0 


2 = 0 


dx2 


UX2 UXl UXl UX2 


where the last equality used the relation in (14.6p . The result now follows since this 
final expression equals det (J(x*)) and since = 1 we have ||w|l 2 > 1. □ 

Theorem 4.6. There exist pi and p 2 in (EH) with a simple root x* € such 

that 


Fx2,Rsyiv)>\\J{F) ^\\l 

and ||J(x *)“^||2 > 1. Thus, an eigenvalue of Rsyiv{x 2 ) can be squared more sensitive 
to perturbations than the corresponding root in the absolute sense. 
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Proof. We give an example for which ||ui ||2 > 1 in Lemma 14.51 For some positive 
parameter u and for some n> 2 consider the polynomials 

Pl{xi,X2) = X 1 X 2 + v}l'^X\, P2{xi,X2) = a~^i{xi + X 2 ) + U^''^X2. 

One can verify that x* = (0, 0) is a common roo1@. Since |6„[(72](0)| = an-ia^\ = 1 
we have ||w ||2 > 1- The result then follows from |det(J(s*))| = and 

Lemma 1331 □ 

Example 4.7. Let us specialize Examvle \3.^ to d = 2, i.e., for some a < 1 and 
q ,2 _|_ ^2 _ 1 consider the system 

Pi = xl + a{axi + (3x2) =0, P 2 = X 2 + a{—j3xi + 0 x 2 ) = 0. 

Again, for the solution = (0,0) we have || Building the Sylvester 

matrix in the monomial basis, we obtain 


RSylv = 


(jj3x2 

X 2 + <Jax2 
0 


(ja 

-a/3 

X 2 + aax2 


1 

0 

-a 13 


As predicted by the theory, X 2 = 0 is an eigenvalue with corresponding right and left 
eigenvectors, respectively, u = [l 0 O] and w = [a/3 aa 1 ] . Moreover, it is 
readily checked that, as expected, R'syivi^)'^ = a^. Therefore, 


k/x 2, Rsylv'j — 


VTT< 


> a 


Theorem 133] mathematically explains the numerical difficulties that practitioners 
have been experiencing with hidden variable resultant methods based on the Sylvester 
resultant. There are successful bivariate rootfinders based on this methodology [39] 
for low degree polynomial systems and it is a testimony to those authors that they 
have developed algorithmic remedies (not cures) for the inherent numerical instability. 

We emphasize that Theorem 14.61 holds for any normalized degree-graded polyno¬ 
mial basis. Thus, the mild numerical instability cannot, in general, be overcome by 
working in a different degree-graded polynomial basis. 

The example in the proof of Theorem 14.61 is quite alarming for a practitioner 
since if u is the unit machine roundoff, then we have || J(0,0)“^||2 = and 

k{x 2, Rsyiv) = u~^. Thus, a numerical rootfinder based on the Sylvester matrix may 
entirely miss a solution that has a condition number larger than A stable 

rootfinder should not miss such a solution. 

When d = 2, we can use Theorem [Q and Lemma 14.51 to conclude that the 
ratio between the conditioning of the Cayley and Sylvester resultant matrices for 
the same eigenvalue X 2 is equal to ||u|| 2 /||rc|| 2 , where v and w are the right and left 
eigenvector of Rsyivix^) associated with the eigenvalue x^- This provides theoretical 
support for the numerical observations in [35] . However, it seems difficult to predict 
a priori if the Cayley or Sylvester matrix will behave better numerically. For real 
polynomials and d = 2, the Cayley resultant matrix is symmetric and this structure 
can be exploited |35j . In the monomial basis, the Sylvester matrix is two stacked 
Toeplitz matrices (see (14.11) 1. It may be that structural differences like these are more 
important than their relatively similar numerical properties when d = 2. 


°By a change of variables, there is an analogous example with a solution anywhere in the complex 
plane. 
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5. A discussion on relative and absolute conditioning. Let X{D) be the 
solution of a mathematical problem depending on data D. In general, with the very 
mild assumption that D and X lie in Banach spaces, it is possible to define the 
absolute condition number of the problem by perturbing the data to -D + SD and 
studying the behaviour of the perturbed solution X{D + 5D) = X{D) + 6X{D, 6D): 


Kabs = lim sup 

<=“^0 ||5D||<e 




Similarly, a relative condition number can be defined by looking at the limit ratios of 
relative changes. 


^rel — lilTL SUp itj:7-\|| Mxaii — ^abs 11 11 • 

'=“^0 ||<5D||<e||D|| ll^-C)|| ||A|| ||A|| 

In this paper, we have compared two absolute condition numbers. One is given 
by Proposition 12.91 there, A = x* is a solution of (HID while D = {pi,... ,pd) is the 
set of polynomials in (HD- The other is given by Lemma 12.111 where is a matrix 
polynomial and X = is the dth component of x*. 

To quote N. J. Higham p. 56]: “Usually, it is the relative condition number 
that is of interest, but it is more convenient to state results for the absolute condition 
number”. This remark applies to our analysis as well. We have found it convenient to 
study the absolute condition number, but when attempting to solve the rootfinding 
problem in floating point arithmetic it is natural to allow for relatively small pertur¬ 
bations, and thus to study the relative condition number. Hence, a natural question is 
whether the exponential increase of the absolute condition number in Theorem l3.7l and 
the squaring in Theorem 14.61 causes a similar effect in the relative condition number. 

It is not immediate that the exponential increase of the absolute condition number 
leads to the same effect in the relative sense. We have found examples where the 
exponential increase of the absolute condition number is perfectly counterbalanced 
by an exponentially small Cayley resultant matrix. For instance, linear polynomial 
systems, when the Cayley resultant method is equivalent to Cramer’s rule, fall into 
this category. In the relative sense, it may be possible to show that the hidden variable 
resultant method based on Cayley or Sylvester is either numerically unstable during 
the construction of the resultant matrix or the resultant matrix has an eigenvalue 
that is more sensitive to small relative perturbations than hoped. We do not know 
yet how to make such a statement precise. 

Instead, we provide an example that shows that the hidden variable resultant 
method remains numerically unstable in the relative sense. Let u be a sufficiently 
small real positive parameter and d > 2. Consider the following polynomial system: 

P 2 i-l{x) = -I- U {^^X 2 i -1 + ^X 2 ?j , 

P2i{x) =xl^+U (^^X2^ - ^X2i-1^ , 1 < f < ld/2\ , 

where if d is odd then take Pd{x) = x"^ + uxd- Selecting H = [—1,1]'^, we have that 
Ibilloo = 1 + for 1 < f < d, except possibly ||pd||oo = 1 -|- u if d is odd. It can be 
shown that the origiiQ, x*, is a simple root, det(J(x*)) = || J(x *)“^||2 = u~^, and 


^By a change of variables, there is an analogous example with a solution anywhere in [—1,1]*^. 
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that 


d-l 

fCayleyi^Sl^ ■ • ■ ? ^d—lj ^1? • ■ • ^ ^d—l') — n (sfe + tk)x\ + (D{u). 

fc=i 

Thus, neither the polynomials pi or the resultant matrix Rcayiey{xd) are small. In 
such an example, the relative condition number will exhibit the same behavior as 
the absolute condition number. In particular, the relative condition number of an 
eigenvalue of Rcayiey(xd) may be larger than the relative condition number of the 
corresponding solution by a factor that grows exponentially with d. 

The same example (for d = 2), and a similar argument, applies to the Sylvester 
matrix showing the conditioning can be squared in the relative sense too. 

6. Future outlook. In this paper we have shown that two popular hidden vari¬ 
able resultant methods based on the Sylvester and Cayley matrices are numerically 
unstable. Our analysis is for degree-graded polynomial bases and does not include 
the Lagrange basis or certain sparse bases. We believe that the analysis of the Cayley 
matrix in Section [3] could be extended to include general polynomial bases, though 
the analysis in Section 3] for the Sylvester matrix is more intimately connected to 
degree-graded bases. We hesitantly suggest that hidden variable resultant methods 
are inherently plagued by numerial instabilities, and that neither other polynomial 
bases nor other resultants can avoid a worst-case scenario that we have identified 
in this paper. We do not know exactly how to formulate such a general statement, 
but we note that practitioners are widely experiencing problems with hidden variable 
resultant methods. In particular, we do not know of a numerical multidimensional 
rootfinder based on resultants that is robust for polynomial systems of large degree n 
and high d. 

However, at the moment the analysis that we offer here is limited to the Cayley 
and Sylvester matrices. Despite our doubts that it exists, we would celebrate the 
discovery of a resultant matrix that can be constructed numerically and that prov- 
ably does not lead to a numerically unstable hidden variable resultant method. This 
would be a breakthrough in global rootfinding with significant practical applications 
as it might allow (HH) to be converted into a large eigenproblem without confronting 
conditioning issues. Solving high-dimensional and large degree polynomial systems 
would then be restricted by computational cost rather than numerical accuracy. 

Finally, we express again our hope that this paper, while appearing rather nega¬ 
tive, will have a positive long-term impact on future research into numerical rootfind- 
ers. 


Acknowledgments. We thank Yuji Nakatsukasa, one of our closest colleagues, 
for his insightful discussions during the writing of [35] that ultimately lead us to con¬ 
sider conditioning issues more closely. We also thank Anthony Austin and Martin 
Lotz for carefully reading a draft and providing us with excellent comments. While 
this manuscript was in a much earlier form Martin Lotz personally sent it to Gre¬ 
gorio Malajovich for his comments. Gregorio’s comprehensive and enthusiastic reply 
encouraged us to proceed with renewed vigor. 

Appendix A. A generalization of Clenshaw’s algorithm for degree- 
graded polynomial bases. This appendix contains the tedious, though necessary, 
proofs required in Section 14.11 for Glenshaw’s algorithm for evaluating polynomials 
expressed in a degree-graded basis. 
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Proof. [Proof of Lemma 14.2] By rearranging (14.31) we have = &fc[p](x) — {akX- 


I3k)bk+i[p]{x) - ELfc+i lj,k+ibj+i[p\ix). Thus, 




p{x) = ao(j)o{x) + '^ 


k^l 


n—1 


bk[p]{x) - {akX + I3k)bk+i[p]{x) - ^ jj^k+ibj+i[p]{x) 

j=k+l 


<j)k{x). 


Now, by interchanging the summations and collecting terms we have 

n n 

p{x) = ao0o(^) + ^ (l>k{x)h[p]{x) - ^ {ak-ix I3k-i) <l>k-i{x)bk[p](x) 


k^l 




n—1 

-E 

i=2 


= aoMx) + (l)i{x)bi[p]{x) 


+ E 

i=i L 

Finally, using we obtain 


b-i 

^ ' 'lj,k+l4>k{x) 

.k=l 


i-1 


(t)j+i{,x) - {ajX + I3j)(t>j{x) - ^-ij^k+i(t>k{x) 


/c=l 


bj+i[pKx) 


bj+i[p]{x) 


n—1 

p{x) = ao(j)o[x) + (j)i{x)bi[p]{x) + ^-ij.i(j)o{x)bj+i[p]{x), 

1=1 

as required. □ 

Section im also shows that Clenshaw’s algorithm connects to the quotient {p{x) — 
p{y))/{x — y). To achieve this we need an immediate result that proves a different re¬ 
currence relation on the Clenshaw shifts to (14.31) . The proof involves tedious algebraic 
manipulations and mathematical strong induction. 

Lemma A.l. Let n be an integer, 0o,..., </>» a degree-graded basis satisfying 621, 
and bn+i[p],... ,bi[p] the Clenshaw shifts satisfying (14.3p . Then, for 1 < j < n, 

n 

bj[(j)n+l]{x) = (UnX I3n)bj[4>n\{x) + ^ ')n,sbj[(j)s-l\{.x). 

s=j+l 


Proof We proceed by induction on j. Let j = n. We have, by (14.31) . 

[0n-t-l] (^) “b /^n)^n+l [^n-t-l] (^) “b [^n] (^); 

where the last equality follows because bn+i[(f>n+i]ix) = bn[(l>n]{x) = 1. Now, suppose 
the result holds for j = n,n — 1,... ,k 1. We have, by (14.31) and the inductive 
hypothesis, 

n 

bk[4>n+i]{x) = {akX-i-I3k)bk+i[(l)n-i-i]{x) + ^ 7j,fe+i6j+i[0„+i](a:) 


j=k-\-l 


= {akX-\- (3k) 


{anX-i-firi)bk-i-l[(l>n]{x) + ^ '^ri,sbk+l[4>s-l]{x) 


s—k-\-2 


n—1 

■ E Tj.fe-Hi 

j=k+l 


(anX/3n)bj+i[(l)n]ix) + ^ 'yn,sbj+l[(l)s-l]ix) 

s^j+2 

'yn,k-\-lbn-\-l [071+1] (^) ■ 
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By interchanging the summations and collecting terms we have 


bk[4>n+lKx) = (anX + Pn) 


s—k-\-3 


{akX + Pk)bk+l[<l)n]ix) + ^ Jj^k+lbj+l[(j)n]ix) 


j=k+l 


{akX + Pk)bk+i[Ps-i]ix) + ^ ■jj^k+ibj+ilPs-iKx) 

Tri,/c+2 [^fc+l] (^) “t” T'n.,/c+l^n+l [0n+l] (^) 


— “t“ Pn')^k[4^n\ H“ ^ ^ ^n,sbk[^s—l\^ 


s—k-\-l 


where in the last equality we used (|4.3D . {akX H- j3k)bk-\-i[4>k-\-i]{x) = bk[(l)k-\-i]{x)^ and 

bn+llpn+lKx) = bk[pk]{x) = 1. □ _ 

The recurrence from Lemma lA.ll allows us to prove Theorem 14.31 
Proof. [Proof of Theorem 14.3| 

Case 1: x y. Since for a fixed y the Clenshaw shifts are linear, i.e., bj[ci(j)i + 
C 24 >k]iy) = cibj[(j)i]{y) + C 2 bj[(j)k]iy) for constants ci and C 2 , it is sufficient to prove 
the theorem for p = (f>n for n > 1. 

We proceed by induction on n. For n = 1 we have 


n—1 

'^ajbj+i[(j)n+i]{y)(l)j = aobi[(j>i]{y) = ao 
3=0 


(fijx) - (fijy) 
x-y 


Assume that the result holds for n = 1,... ,k — 1. From the inductive hypothesis, 
we have 


Pk+ijx) - pk+ijy) 
x-y 


oikPki.x) + {akX + 

x-y 

+ - 

k-l 

akPkix) + {akX + Pk) ctjbj+i [Pk] iy)4>j jx) 

3=0 


k 3-2 

+ '^lk,j'^asbs+i[(t)j-i]{y)(j)s{x). 
j=l s=0 


Moreover, by interchanging the summations and collecting terms we have 


Pk+lix) (fk+liy) 
X-y 


= akpkix) + {akX + Pk)ak-ibk[<l)k]iy)4>k-i{x) 


k-2 

+ E' 

3=0 


{akX + Pk)bj+i[(j)k\{y) + ^ -^k,sbj+i[(j)s-i]{y) 

s=3+2 


pj{x). 


Finally, since bk+i[4>k+i]{y) = 1, 6fe[<^fe+i](y) = [ukX + Pk)bk[4>k\{y), and by ((Ol) . we 
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have 


^ + ak-ihk[(l)k+i\{y)4)k-i{x) 

X y 

k-2 

+ J2 ^j^j+^i^k+i]{y)(l)j{x) 

J=0 


and the result follows by induction. 

Case 2: x = y. Immediately follows from x y hy using L’Hospital’s rule 
on (|4.5I) . □ 


REFERENCES 

[1] E. L. Allgower, K. Georg, and R. Miranda, The method of resultants for computing real 

solutions of polynomial systems, SIAM J. Numer. Anal., 29 (1992), pp. 831-844. 

[2] J. Asakura, T. Sakurai, H. Tadano, T. Ikegami, and K. Kimura, a numerical method for 

polynomial eigenvalue problems using contour integral, Japan J. Indust. Appl. Math., 27 
(2010), pp. 73-90. 

[3] D. J. Bates, J. D. Hauenstein, A. J. Sommese, and C. W. Wampler, Numerically Solving 

Polynomial Systems with Bertini, SIAM, 2013. 

[4] D. A. Bini and a. Marco, Computing curve intersection by means of simultaneous iterations, 

Numer. Algor., 43 (2006), pp. 151-175. 

[5] D. A. Bini and V. Noferini, Solving polynomial eigenvalue problems by means of the Ehrlich- 

Aberth method. Linear Algebra Appl., 439 (2013), pp. 1130-1149. 

[6] D. Bondyfalat, B. Mourrain, and V. Y. Pan, Solution of a polynomial system of equations 

via the eigenvector computation, Linear Algebra Appl., 319 (2000), pp. 193-209. 

[7] J. P. Boyd, Computing zeros on a real interval through Chebyshev expansion and polynomial 

rootfinding, SIAM J. Numer. Anal., 40 (2002), pp. 1666-1682. 

[8] J. P. Boyd, Solving Transcendental Equations: The Chebyshev Polynomial Proxy and Other 

Numerical Rootfinders, SIAM, 2014. 

[9] B. BuCHBERGER, An algorithm for finding the basis elements of the residue class ring of a zero 

dimensional polynomial ideal, J. Symbolic Comput., 41 (2006), pp. 475-511. 

[10] L. Buse, H. Khalil, and B. Mourrain, Resultant-based methods for plane curves intersection 

problems. Computer algebra in scientific computing. Springer Berlin Heidelberg, 3718 (2005), 
pp. 75-92. 

[11] E. Cattani and a. Dickenstein, Introduction to residues and resultants, in, A. Dickenstein 

AND 1. Z. Emiris, editors. Solving Polynomial Equations. Foundations, Algorithms, and 
Applications, Springer, 2005. 

[12] A. Cayley, On the theory of elimination, Cambridge and Dublin Math. J. Ill, (1848), pp. 116— 

120 . 

[13] E.-W. Chionh, M. Zhang, and R. N. Goldman, Fast computation of the Bezout and Dixon 

resultant matrices, J. Symb. Comput., 33 (2002), pp. 13-29. 

[14] C. W. Clenshaw, a note on the summation of Chebyshev series. Math. Comput., 9 (1955), 

pp. 118-120. 

[15] D. A. Cox, J. Little, and D. O’Shea, Using Algebraic Geometry, Springer, 2013. 

[16] F. De Teran, F. M. Dopico, and D. S. Mackey, Fiedler companion linearizations and the 

recovery of minimal indices, SIAM J. Mat. Anal. Appl., 31 (2010), pp. 2181-2204. 

[17] P. Dreesen, K. Batselier, and B. De Moor, Back to the roots: Polynomial system solv¬ 

ing, linear algebra, systems theory, Proc. 16th IFAC Symposium on System Identification 
(SYSID), 2012, pp. 1203-1208. 

[18] 1. Z. Emiris, Sparse elimination and applications in kinematics. Dissertation, University of 

California, Berkeley, 1994. 

[19] 1. Z. Emiris and V. Z. Pan, Symbolic and numeric methods for exploiting structure in con¬ 

structing resultant matrices, J. Symbolic Comput., 33 (2002), pp. 393—413. 

[20] 1. M. Gelfand, M. Kapranov, and A. Zelevinsky, Discriminants, Resultants, and Multidi¬ 

mensional Determinants, Springer, Birkhauser, Boston, 2008. 

[21] L. Gemignani, and V. Noferini, The Ehrlich-Aberth method for palindromic matrix polyno¬ 

mials represented in the Dickson basis. Linear Algebra Appl., 438 (2013), pp. 1645—1666. 



24 


[22] I. Gohberg, P. Lancaster, and L. Rodman, Matrix Polynomials, SIAM, Philadelphia, USA, 

2009, (unabridged republication of book first published by Academic Press in 1982). 

[23] I. J. Good, The colleague matrix, a Chebyshev analogue of the companion matrix, The Quar¬ 

terly Journal of Mathematics, 12 (1961), pp. 61-68. 

[24] N. J. Hicham, Accuracy and Stability of Numerical Algorithms, SIAM, 2nd edition, 2002. 

[25] N. J. Hicham, Function of Matrices: Theory and Computation, SIAM, 2008. 

[26] R. A. Horn, and C. R. Johnson, Matrix Analysis, Cambridge University Press, 2nd edition, 

New York, 2013. 

[27] G. JONSSON AND S. VavasiS, Accurate solution of polynomial equations using Macaulay resultant 

matrices, Math. Comput., 74 (2005), pp. 221-262. 

[28] D. Kapur and T. Saxena, Comparison of various multivariate resultant formulations, ACM 

Proceedings of the 1995 international symposium on Symbolic and algebraic computation, 
1995. 

[29] F. C. Kirwan, Complex Algebraic Curves, Cambridge University Press, Cambridge, 1992. 

[30] H. Li, A simple solution to the six-point two-view focal-length problem, Computer VisionECCV, 

Springer Berlin Heidelberg, (2006), pp. 200-213. 

[31] D. S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann, Vector spaces of linearizations for 

matrix polynomials, SIAM J. Matrix Anal. AppL, 28 (2006), pp. 971-1004. 

[32] D. Manocha and J. F. Canny, Multipolynomial resultants and linear algebra. Papers from the 

international symposium on Symbolic and algebraic computation. ACM, 1992. 

[33] D. Manocha and J. Demmel, Algorithms for intersecting parametric and algebraic curves I: 

simple intersections, ACM Trans. Graphics, 13 (1994), pp. 73-100. 

[34] Y. Nakatsukasa, V. Noferini, and A. Townsend, Vector spaces of linearizations for matrix 

polynomials: a bivariate polynomial approach, submitted. 

[35] Y. Nakatsukasa, V. Noferini, and A. Townsend, Computing the common zeros of two 

bivariate functions via Bezout resultants, Numer. Math., 129 (2015), pp. 181—209. 

[36] S. Ragnarsson and C. F. Van Loan, Block tensor unfoldings, SIAM J. Mat. Anal. Apph, 33 

(2012), pp. 149-169. 

[37] M.-R. Skrzipek, Polynomial evaluation and associated polynomials, Numer. Math., 79 (1998), 

pp. 601-613. 

[38] A. J. SOMMESE AND C. W. Wampler, The Numerical Solution of Systems of Polynomials 

Arising in Engineering and Science, World Scientific, 2005. 

[39] L. SORBER, M. Van Barel, and L. De Lathauwer, Numerical solution of bivariate and 

polyanalytic polynomial systems, SIAM J. Numer. Anal., 52 (2014), pp. 1551—1572. 

[40] M. 1. Syam, Finding all real zeros of polynomial systems using multi-resultant, J. Comput. 

Appl. Math., 167 (2004), pp. 417-428. 

[41] L. Taslaman, An algorithm for quadratic eigenproblems with low rank damping , SIAM J. 

Matrix Anal. Appl., 36 (2015), pp. 251—272. 

[42] F. Tisseur, Backward error and condition of polynomial eigenvalue problems. Linear Algebra 

Appl., 309 (2000), pp. 339-361. 

[43] F. Tisseur and K. Meerbergen, The quadratic eigenvalue problem, SIAM Review, 43 (2001), 

pp. 235-286. 

[44] A. Townsend, Computing with functions of two variables, DPhil Thesis, The University of 

Oxford, 2014. 

[45] A. Townsend and L. N. Trefethen, An extension of Chebfun to two dimensions, SIAM J. 

Sci. Comput., 35 (2013), C495-C518. 

[46] J. Weiss, Resultant methods for the inverse kinematics problem. Computational kinematics. 

Springer Netherlands, 28 (1993), pp. 41-52. 



